#include<../data_tools/usefull_routines.c>

size_t *ind;
double *height;
double *radio;
particulas *diskpart;
int ndisk;

int DensityVertical(int N_bines_vertical, int imin, int imax, int Npuntos, FILE *parametros);


int main(int argc, char *argv[])
{
  int i,j,nreal,nbines,lim[2],N_total,npart,status,nbinesv,limpart;
  float rmin,rmax,deltar,rini,rfin,area;
  float mass,meanmass;
  double *R,*SigmaR,*errors,f,p,size;
  char *infile,outfiles[1000];
 
  infile = argv[1];  
  nbines = atoi(argv[2]);  
  nbinesv = atoi(argv[3]);
  limpart = 100;  
  //limpart = atoi(argv[4]);  

  printf("\n");
  printf("#############################################################\n");
  printf("Computing density profile for the disk:\n");
  printf("%s\n",argv[1]);
  printf("#############################################################\n");
  printf("\n");

  //Reading datas
  FILE *diskfile;
  diskfile = fopen(infile,"r");
  ndisk = contadorParticulas(infile);
  printf("Reading %d particles\n",ndisk);
  
  diskpart = (particulas *)malloc((size_t)ndisk*sizeof(particulas));
  if(diskpart == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  radio = (double *)malloc((size_t)ndisk*sizeof(double));
  if(radio == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  height = (double *)malloc((size_t)ndisk*sizeof(double));
  if(height == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  R = (double *)malloc((size_t)ndisk*sizeof(double));
  if(R == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  SigmaR = (double *)malloc((size_t)ndisk*sizeof(double));
  if(SigmaR == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  errors = (double *)malloc((size_t)ndisk*sizeof(double));
  if(errors == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

 ind = (size_t *)malloc((size_t)ndisk*sizeof(size_t));
  if(ind == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  for(i=0; i<ndisk; i++)
    {
      fscanf(diskfile,"%d %e %e %e %e %e %e %e %e\n",
	     &diskpart[i].id,
	     &diskpart[i].pos[0],&diskpart[i].pos[1],&diskpart[i].pos[2],
	     &diskpart[i].vel[0],&diskpart[i].vel[1],&diskpart[i].vel[2],
	     &diskpart[i].mass,
	     &diskpart[i].pot);
      radio[i] = sqrt(pow(diskpart[i].pos[0],2) + pow(diskpart[i].pos[1],2));
      height[i] = diskpart[i].pos[2];
      ind[i] = i;
    }
 
  gsl_sort_index(ind,radio,1,ndisk);

  rmin = log10(radio[ind[0]]);
  //rmax = log10(radio[ind[ndisk-1]]);
  rmax = log10(10.0);
  printf("Done.\n");
  printf("\n");


  //Computing radial density profile
  printf("Computing radial density profile\n");
  FILE *perfilfile;
  sprintf(outfiles,"%s_radial_density.dat",infile);
  perfilfile = fopen(outfiles,"w");
 
  deltar = (rmax-rmin)/nbines;
  rini = rmin;
  N_total = 0;
  nreal = 0;
  for(i=0; i<nbines; i++)
    {
      rfin = rini+deltar;
      npart = 0;
      mass = 0.0;
      lim[0] = lim[1] = 0;

      for(j=0; j<ndisk; j++)
	{
	  if((radio[ind[j]]>=pow(10,rini)) && (radio[ind[j]]<pow(10,rfin)))
	    {
	      npart = npart+1;
	      if(npart==1)
		{
		  lim[0]=j;
		}
	      mass = mass+diskpart[ind[j]].mass;
	      lim[1]=j;
	    }
	}
      

      if(npart>=limpart)
	{
	  area=2.0*M_PI*pow(10,rini)*(pow(10,rfin)-pow(10,rini));
	  meanmass = mass/npart;
	  R[i] = pow(10,rini);
	  SigmaR[i] = mass/area; 
	  errors[i] = meanmass*sqrt(npart)/area;
	  fprintf(perfilfile,"%e %e %e\t",R[i],SigmaR[i],errors[i]);
	  nreal = nreal+1;
	  DensityVertical(nbinesv,lim[0],lim[1],npart,perfilfile);
	}
      if(npart<limpart)
	{
	  R[i] = 0.0;
	  SigmaR[i] = 0.0; 
	  errors[i] = 1.0;
	}

		
      //seguir acá
      N_total = N_total+npart;
      rini = rfin;
    }
  printf("Done.\n");
  printf("\n");


 //Fitting density profile
  printf("Fitting density profile\n");
  double densidadRadial(double R,double h,double Md)
  {
    double perfil;
    
    perfil = Md*(1.0/(4.0*M_PI*h*h))*exp(-R/h);

    return perfil;
  }
  //funcion que calcula el chisquare
  double ChiSquare(const gsl_vector *parmodelo, void *param)
  {
    int i;
    double h,Md;
    double chisquare;
    double xo,yo,s,yt;
    
    h = gsl_vector_get(parmodelo,0);
    Md = gsl_vector_get(parmodelo,1);
    
    chisquare=0;
    for(i=0; i<nbines; i++)
      {
	if(R[i]!=0.0)
	  {
	    xo=R[i];
	    yo=SigmaR[i];
	    s=errors[i];
	    yt=densidadRadial(xo,h,Md);
	    chisquare+=gsl_pow_int(yt-yo,2)/(s*s);
	  }
      }
    
    return chisquare;
  }


  //--------------------------
  //MINIMIZACION DEL CHISQUARE
  gsl_vector *parmodelo=gsl_vector_alloc(2);
  gsl_vector *dpar=gsl_vector_alloc(2);
  
  gsl_vector_set(parmodelo,0,1.0);
  gsl_vector_set(parmodelo,1,1.0);
  
  gsl_vector_set_all(dpar,0.1);
  
  //1-DECLARACION
  const gsl_multimin_fminimizer_type *tipo;
  gsl_multimin_fminimizer *minimizador;
  gsl_multimin_function funcion;
  
  //2-ALOCACION E INICIALIZACION
  tipo=gsl_multimin_fminimizer_nmsimplex;
  minimizador=
    gsl_multimin_fminimizer_alloc(tipo,2);
  
  funcion.f=&ChiSquare;
  funcion.n=2;
  funcion.params=NULL;
  
  gsl_multimin_fminimizer_set(minimizador,
			      &funcion,
			      parmodelo,
			      dpar
			      );
  
  //3-USO
  do{
    gsl_multimin_fminimizer_iterate(minimizador);
    //Es un estimativo de la distancia al minimo
    size=gsl_multimin_fminimizer_size(minimizador);
    status=gsl_multimin_test_size(size,1E-5);
  }while(status!=GSL_SUCCESS);
  
  parmodelo=gsl_multimin_fminimizer_x(minimizador);
  fprintf(stdout,"h = %e\n",
	  gsl_vector_get(parmodelo,0));
  fprintf(stdout,"Md = %e\n",
	  gsl_vector_get(parmodelo,1));  
  fprintf(stdout,"Chisquare = %e with %d points\n",
	  gsl_multimin_fminimizer_minimum(minimizador),nreal);
  printf("Done.\n");
  printf("\n");

  FILE *fitfile;
  sprintf(outfiles,"%s_radial_density.fit",infile);
  fitfile = fopen(outfiles,"w");
  
  for(f=pow(10,rmin); f<pow(10,rmax); f+=0.01)
    {
      p = densidadRadial(f,gsl_vector_get(parmodelo,0),gsl_vector_get(parmodelo,1));
      fprintf(fitfile,"%e %e\n",f,p);
      
    }
  
   
  //4-LIBERACION
  gsl_multimin_fminimizer_free(minimizador);
  

  free(diskpart);
  free(R);
  free(SigmaR);
  free(errors);
  free(ind);

  fclose(perfilfile);
  fclose(diskfile);
  fclose(fitfile);

  return 0;
}


int DensityVertical(int N_bines_vertical, int imin, int imax, int Npuntos, FILE *parametros)
{
  int i,j,status,N_total,npart,N_tomados;
  size_t *index;  
  
  double *Z,*Sigma,*errores,*densidad_numero,*z;
  double zmin,zmax,delta_z,z_ini,z_fin,masa,linea;  
  double size,masamed; 

  //calcula la densidad
  //carga los datos
  Z=malloc(Npuntos*sizeof(double));
  z=malloc(N_bines_vertical*sizeof(double));
  index=malloc(Npuntos*sizeof(size_t));
  Sigma=malloc(N_bines_vertical*sizeof(double));
  densidad_numero=malloc(N_bines_vertical*sizeof(double));
  errores=malloc(N_bines_vertical*sizeof(double));

  //printf("Numero de puntos %d %d\n",Npuntos,sizeof(Z));
  i=0;
  N_tomados=N_bines_vertical;
  for(j=0;j<ndisk;j++) 
    {
      if((j>=imin) && (j<=imax))
	{
	  Z[i] = height[ind[j]];
	  index[i] = ind[j];
	  i=i+1;
	}
    }
  
  //calcula el perfil de densidad vertical del bin
  gsl_sort_index(index,Z,1,Npuntos);

  zmin=(-1.0)*Z[index[Npuntos-1]];
  zmax=Z[index[Npuntos-1]];

  delta_z=(zmax-zmin)/N_bines_vertical; // mirar como funciona
  z_ini=zmin;

  N_total=0;

  for(i=0; i<N_bines_vertical; i++)
    {
      z_fin=z_ini+delta_z;
      npart=0;
      masa=0.0;
      for(j=0; j<Npuntos; j++)
	{
	  if((Z[j]>=z_ini) && (Z[j]<z_fin))
	    {
	      npart=npart+1;
	      masa=masa+diskpart[index[j]].mass;
	    }
	}
      N_total=N_total+npart;

      linea=z_fin-z_ini;

      masamed=masa/npart;
      Sigma[i]=masa/linea;
      errores[i]=masamed*sqrt(npart)/linea;
      densidad_numero[i]=npart/linea;
      z[i]=z_ini;   

      if(npart==0) 
	{
	  Sigma[i]=0.0;
	  z[i]=0.0;
	  errores[i]=1.0;
	  N_tomados=N_tomados-1;;
	}

      z_ini=z_fin;
    }
  printf("numero total de particulas encontradas en el bin %d\n",N_total);
  printf("N_tomados=%d\n",N_tomados);

  //ajusta el perfil calculado arriba minimizando el chisquare

  //funcion esperada para el perfil
  double densidadVertical(double z,double sigma_r,double zo,double zmax)
  {
    double Cosh,perfil;
    
    Cosh = 1.0/(cosh((z-zmax)/zo));
    
    perfil = (sigma_r)*(1.0/zo)*(pow(Cosh,2));
    
    return perfil;
  }
  //funcion que calcula el chisquare
  double ChiSquare(const gsl_vector *parmodelo, void *param)
  {
    int i;
    double sigma_r,zmax,zo;
    double chisquare;
    double xo,yo,s,yt;
    
    sigma_r=gsl_vector_get(parmodelo,0);
    zo=gsl_vector_get(parmodelo,1);
    zmax=gsl_vector_get(parmodelo,2);
    
    //CALCULO DEL CHI-SQUARE
    chisquare=0;
    for(i=0;i<N_bines_vertical;i++)
      {
	//Calculo valores teoricos
	xo=z[i];
	yo=Sigma[i];
	s=errores[i];
	yt=densidadVertical(xo,sigma_r,zo,zmax);
	//Calculo del chisquare
	chisquare+=gsl_pow_int(yt-yo,2)/(s*s);
      }
    
    return chisquare;
  }
  //--------------------------
  //MINIMIZACION DEL CHISQUARE
  gsl_vector *parmodelo=gsl_vector_alloc(3);
  gsl_vector *dpar=gsl_vector_alloc(3);
  
  gsl_vector_set(parmodelo,0,1.0);
  gsl_vector_set(parmodelo,1,1.0);
  gsl_vector_set(parmodelo,2,0.0);
  
  gsl_vector_set_all(dpar,0.1);
  
  //1-DECLARACION
  const gsl_multimin_fminimizer_type *tipo;
  gsl_multimin_fminimizer *minimizador;
  gsl_multimin_function funcion;
  
  //2-ALOCACION E INICIALIZACION
  tipo=gsl_multimin_fminimizer_nmsimplex;
  minimizador=
    gsl_multimin_fminimizer_alloc(tipo,3);
  
  funcion.f=&ChiSquare;
  funcion.n=3;
  funcion.params=NULL;
  
  gsl_multimin_fminimizer_set(minimizador,
			      &funcion,
			      parmodelo,
			      dpar
			      );
  
  //3-USO
  do{
    gsl_multimin_fminimizer_iterate(minimizador);
    //Es un estimativo de la distancia al minimo
    size=gsl_multimin_fminimizer_size(minimizador);
    status=gsl_multimin_test_size(size,1E-5);
  }while(status!=GSL_SUCCESS);
  
  parmodelo=gsl_multimin_fminimizer_x(minimizador);
  /*fprintf(stdout,"masa_r = %e\n",
	  gsl_vector_get(parmodelo,0));
  fprintf(stdout,"zo = %e\n",
	  gsl_vector_get(parmodelo,1));
  fprintf(stdout,"zmax = %e\n",
	  gsl_vector_get(parmodelo,2));
  fprintf(stdout,"Chisquare = %e\n",
	  gsl_multimin_fminimizer_minimum(minimizador));
  */
  /*
  for(f=zmin;f<zmax;f+=0.01)
    {
      p=densidadVertical(f,gsl_vector_get(parmodelo,0),gsl_vector_get(parmodelo,1),gsl_vector_get(parmodelo,2));
      //  fprintf(pf,"%e %e \n",f,p);
   
    }*/
 
  fprintf(parametros,"%e\t%e\t%e\n",
	  gsl_vector_get(parmodelo,0),
	  gsl_vector_get(parmodelo,1),
	  gsl_vector_get(parmodelo,2));

  
  //4-LIBERACION
  gsl_multimin_fminimizer_free(minimizador);
  
  free(Z);
  free(z);
  free(index);
  free(Sigma);
  free(densidad_numero);
  free(errores);

  return 0;
  
}
